<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<!-- saved from url=(0014)about:internet -->
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">

<title>Chapter 6: Exercise 9</title>

<base target="_blank"/>

<style type="text/css">
body, td {
   font-family: sans-serif;
   background-color: white;
   font-size: 12px;
   margin: 8px;
}

tt, code, pre {
   font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
}

h1 { 
   font-size:2.2em; 
}

h2 { 
   font-size:1.8em; 
}

h3 { 
   font-size:1.4em; 
}

h4 { 
   font-size:1.0em; 
}

h5 { 
   font-size:0.9em; 
}

h6 { 
   font-size:0.8em; 
}

a:visited {
   color: rgb(50%, 0%, 50%);
}

pre {	
   margin-top: 0;
   max-width: 95%;
   border: 1px solid #ccc;
}

pre code {
   display: block; padding: 0.5em;
}

code.r {
   background-color: #F8F8F8;
}

table, td, th {
  border: none;
}

blockquote {
   color:#666666;
   margin:0;
   padding-left: 1em;
   border-left: 0.5em #EEE solid;
}

hr {
   height: 0px;
   border-bottom: none;
   border-top-width: thin;
   border-top-style: dotted;
   border-top-color: #999999;
}

@media print {
   * { 
      background: transparent !important; 
      color: black !important; 
      filter:none !important; 
      -ms-filter: none !important; 
   }

   body { 
      font-size:12pt; 
      max-width:100%; 
   }
       
   a, a:visited { 
      text-decoration: underline; 
   }

   hr { 
      visibility: hidden;
      page-break-before: always;
   }

   pre, blockquote { 
      padding-right: 1em; 
      page-break-inside: avoid; 
   }

   tr, img { 
      page-break-inside: avoid; 
   }

   img { 
      max-width: 100% !important; 
   }

   @page :left { 
      margin: 15mm 20mm 15mm 10mm; 
   }
     
   @page :right { 
      margin: 15mm 10mm 15mm 20mm; 
   }

   p, h2, h3 { 
      orphans: 3; widows: 3; 
   }

   h2, h3 { 
      page-break-after: avoid; 
   }
}

</style>

<!-- Styles for R syntax highlighter -->
<style type="text/css">
   pre .operator,
   pre .paren {
     color: rgb(104, 118, 135)
   }

   pre .literal {
     color: rgb(88, 72, 246)
   }

   pre .number {
     color: rgb(0, 0, 205);
   }

   pre .comment {
     color: rgb(76, 136, 107);
   }

   pre .keyword {
     color: rgb(0, 0, 255);
   }

   pre .identifier {
     color: rgb(0, 0, 0);
   }

   pre .string {
     color: rgb(3, 106, 7);
   }
</style>

<!-- R syntax highlighter -->
<script type="text/javascript">
var hljs=new function(){function m(p){return p.replace(/&/gm,"&amp;").replace(/</gm,"&lt;")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
hljs.initHighlightingOnLoad();
</script>


<!-- MathJax scripts -->
<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>



</head>

<body>
<h1>Chapter 6: Exercise 9</h1>

<p>Load the Boston dataset</p>

<pre><code class="r">set.seed(1)
library(MASS)
attach(Boston)
</code></pre>

<h3>a</h3>

<pre><code class="r">lm.fit = lm(nox ~ poly(dis, 3), data = Boston)
summary(lm.fit)
</code></pre>

<pre><code>## 
## Call:
## lm(formula = nox ~ poly(dis, 3), data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12113 -0.04062 -0.00974  0.02338  0.19490 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(&gt;|t|)    
## (Intercept)    0.55470    0.00276  201.02  &lt; 2e-16 ***
## poly(dis, 3)1 -2.00310    0.06207  -32.27  &lt; 2e-16 ***
## poly(dis, 3)2  0.85633    0.06207   13.80  &lt; 2e-16 ***
## poly(dis, 3)3 -0.31805    0.06207   -5.12  4.3e-07 ***
## ---
## Signif. codes:  0 &#39;***&#39; 0.001 &#39;**&#39; 0.01 &#39;*&#39; 0.05 &#39;.&#39; 0.1 &#39; &#39; 1 
## 
## Residual standard error: 0.0621 on 502 degrees of freedom
## Multiple R-squared: 0.715,   Adjusted R-squared: 0.713 
## F-statistic:  419 on 3 and 502 DF,  p-value: &lt;2e-16
</code></pre>

<pre><code class="r">dislim = range(dis)
dis.grid = seq(from = dislim[1], to = dislim[2], by = 0.1)
lm.pred = predict(lm.fit, list(dis = dis.grid))
plot(nox ~ dis, data = Boston, col = &quot;darkgrey&quot;)
lines(dis.grid, lm.pred, col = &quot;red&quot;, lwd = 2)
</code></pre>

<p><img src="" alt="plot of chunk 9a"/> </p>

<p>Summary shows that all polynomial terms are significant while predicting nox using dis. Plot shows a smooth curve fitting the data fairly well.</p>

<h3>b</h3>

<p>We plot polynomials of degrees 1 to 10 and save train RSS.</p>

<pre><code class="r">all.rss = rep(NA, 10)
for (i in 1:10) {
    lm.fit = lm(nox ~ poly(dis, i), data = Boston)
    all.rss[i] = sum(lm.fit$residuals^2)
}
all.rss
</code></pre>

<pre><code>##  [1] 2.769 2.035 1.934 1.933 1.915 1.878 1.849 1.836 1.833 1.832
</code></pre>

<p>As expected, train RSS monotonically decreases with degree of polynomial. </p>

<h3>c</h3>

<p>We use a 10-fold cross validation to pick the best polynomial degree.</p>

<pre><code class="r">library(boot)
all.deltas = rep(NA, 10)
for (i in 1:10) {
    glm.fit = glm(nox ~ poly(dis, i), data = Boston)
    all.deltas[i] = cv.glm(Boston, glm.fit, K = 10)$delta[2]
}
plot(1:10, all.deltas, xlab = &quot;Degree&quot;, ylab = &quot;CV error&quot;, type = &quot;l&quot;, pch = 20, 
    lwd = 2)
</code></pre>

<p><img src="" alt="plot of chunk 9c"/> </p>

<p>A 10-fold CV shows that the CV error reduces as we increase degree from 1 to 3, stay almost constant till degree 5, and the starts increasing for higher degrees. We pick 4 as the best polynomial degree.</p>

<h3>d</h3>

<p>We see that dis has limits of about 1 and 13 respectively. We split this range in roughly equal 4 intervals and establish knots at \( [4, 7, 11] \). Note: bs function in R expects either df or knots argument. If both are specified, knots are ignored.</p>

<pre><code class="r">library(splines)
sp.fit = lm(nox ~ bs(dis, df = 4, knots = c(4, 7, 11)), data = Boston)
summary(sp.fit)
</code></pre>

<pre><code>## 
## Call:
## lm(formula = nox ~ bs(dis, df = 4, knots = c(4, 7, 11)), data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1246 -0.0403 -0.0087  0.0247  0.1929 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(&gt;|t|)
## (Intercept)                             0.7393     0.0133   55.54  &lt; 2e-16
## bs(dis, df = 4, knots = c(4, 7, 11))1  -0.0886     0.0250   -3.54  0.00044
## bs(dis, df = 4, knots = c(4, 7, 11))2  -0.3134     0.0168  -18.66  &lt; 2e-16
## bs(dis, df = 4, knots = c(4, 7, 11))3  -0.2662     0.0315   -8.46  3.0e-16
## bs(dis, df = 4, knots = c(4, 7, 11))4  -0.3980     0.0465   -8.56  &lt; 2e-16
## bs(dis, df = 4, knots = c(4, 7, 11))5  -0.2568     0.0900   -2.85  0.00451
## bs(dis, df = 4, knots = c(4, 7, 11))6  -0.3293     0.0633   -5.20  2.9e-07
##                                          
## (Intercept)                           ***
## bs(dis, df = 4, knots = c(4, 7, 11))1 ***
## bs(dis, df = 4, knots = c(4, 7, 11))2 ***
## bs(dis, df = 4, knots = c(4, 7, 11))3 ***
## bs(dis, df = 4, knots = c(4, 7, 11))4 ***
## bs(dis, df = 4, knots = c(4, 7, 11))5 ** 
## bs(dis, df = 4, knots = c(4, 7, 11))6 ***
## ---
## Signif. codes:  0 &#39;***&#39; 0.001 &#39;**&#39; 0.01 &#39;*&#39; 0.05 &#39;.&#39; 0.1 &#39; &#39; 1 
## 
## Residual standard error: 0.0619 on 499 degrees of freedom
## Multiple R-squared: 0.718,   Adjusted R-squared: 0.715 
## F-statistic:  212 on 6 and 499 DF,  p-value: &lt;2e-16
</code></pre>

<pre><code class="r">sp.pred = predict(sp.fit, list(dis = dis.grid))
plot(nox ~ dis, data = Boston, col = &quot;darkgrey&quot;)
lines(dis.grid, sp.pred, col = &quot;red&quot;, lwd = 2)
</code></pre>

<p><img src="" alt="plot of chunk 9d"/> </p>

<p>The summary shows that all terms in spline fit are significant. Plot shows that the spline fits data well except at the extreme values of \( dis \), (especially \( dis > 10 \)). </p>

<h3>e</h3>

<p>We fit regression splines with dfs between 3 and 16. </p>

<pre><code class="r">all.cv = rep(NA, 16)
for (i in 3:16) {
    lm.fit = lm(nox ~ bs(dis, df = i), data = Boston)
    all.cv[i] = sum(lm.fit$residuals^2)
}
all.cv[-c(1, 2)]
</code></pre>

<pre><code>##  [1] 1.934 1.923 1.840 1.834 1.830 1.817 1.826 1.793 1.797 1.789 1.782
## [12] 1.782 1.783 1.784
</code></pre>

<p>Train RSS monotonically decreases till df=14 and then slightly increases for df=15 and df=16.</p>

<h3>f</h3>

<p>Finally, we use a 10-fold cross validation to find best df. We try all integer values of df between 3 and 16.</p>

<pre><code class="r">all.cv = rep(NA, 16)
for (i in 3:16) {
    lm.fit = glm(nox ~ bs(dis, df = i), data = Boston)
    all.cv[i] = cv.glm(Boston, lm.fit, K = 10)$delta[2]
}
</code></pre>

<pre><code>## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
## Warning: some &#39;x&#39; values beyond boundary knots may cause ill-conditioned bases
</code></pre>

<pre><code class="r">plot(3:16, all.cv[-c(1, 2)], lwd = 2, type = &quot;l&quot;, xlab = &quot;df&quot;, ylab = &quot;CV error&quot;)
</code></pre>

<p><img src="" alt="plot of chunk 9f"/> </p>

<p>CV error is more jumpy in this case, but attains minimum at df=10. We pick \( 10 \) as the optimal degrees of freedom.</p>

</body>

</html>

